#include "Tria3ShellCalculator.h"
#include <iostream>

#include "Utility/Logging/SGLogger.h"
#include "DataStructure/Input/Solver/DOFDirections.h"

#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Result/Nodal/IDISPLACEMENTService.h"
#include "DBManager/IService/Result/IntegStress/ICTRIA3IntegStressService.h"

#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/CoordinateProcessor.h"
#include "SGFem/ElementCalculator/Common/ElementTools.hpp"
#include "SGFem/ShellSecCalculator/ShellSectionParser1.h"

#include "Operater.h"

using namespace std;
using namespace SG::Algebra;
using SG::DataStructure::Common::PI;
using namespace SG::DataStructure::Common;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::CoordinateProcessor;
using namespace SG::FEMSolver::Element;
using SG::FEMSolver::Element::ShellSectionParser1;

Tria3ShellCalculator::Tria3ShellCalculator()
{
}

void Tria3ShellCalculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type = type;
    m_mem.Initialize ();
    m_plate.Initialize();
}

/** @brief 单元刚度矩阵计算
 * 
 * 分别调用膜单元和板单元计算器对刚度矩阵贡献。按照位置进行组装后将其旋转到全局坐标系下
 */
void Tria3ShellCalculator::ComputeStiffness (const Matrix& coord, const ShellSectionProp1& secConstant, const Real zoff, const Matrix& Tele, Matrix& _OUT kg)
{
    static const std::vector<std::size_t> indexm {0, 1, 5, 6, 7, 11, 12, 13, 17};
    /** ToDo 变厚度壳需要对中心点进行插值 */
    auto tp = secConstant.m_t;
    Real thickness[3] = {tp, tp, tp};

    Matrix keLocal (18, 18);
    Matrix Bm (3, 9);
    if (false == secConstant.nomembrane)
    {
        const auto& A = secConstant.m_A;
        Matrix Cm (3,3, {
            A(0,0), A(1,0), A(2,0),
            A(1,0), A(3,0), A(4,0),
            A(2,0), A(4,0), A(5,0)
        });
        Cm = Cm* tp;

        // 调用Membrane单元进行计算
        Matrix km (9, 9);
        m_mem.ComputeStrainMatrix (coord, _OUT Bm);

        const Real& x2 = coord(1, 0);
        const Real& x3 = coord(2, 0);
        const Real& y3 = coord(2, 1);
        const Real A2 = y3 * x2;
        auto BT = Bm.Transpose();
        km = BT * Cm * Bm * (A2 * 0.5);

        Matrix Q (1, 9,{- (x2 - x3) / (A2 * 2.0),  y3 / (A2 * 2.0), 1.0 / 3.0, 
                        -x3 / (A2 * 2.0),       -y3 / (A2 * 2.0), 1.0 / 3.0, 
                        x2 / (A2 * 2.0),                    0.0, 1.0 / 3.0});

        auto GV = A2 / 2.0 * Cm (2, 2);

        // drill
        const Real KROT6 = 100.0;
        Real delta2{1.0e-4 * KROT6};
        Real tmp = delta2 * GV;
        Matrix Sr =  Q.Transpose() * (Q * tmp);

        km += Sr;
        auto lambda =  1.0e-10 * GV * KROT6;
        km(2, 2) += lambda;
        km(5, 5) += lambda;
        km(8, 8) += lambda;

        keLocal.BlockFill (indexm, indexm, km);
    }

    Matrix kp (9, 9);
    Matrix membPlateKe (9, 9);          //拉弯刚度
    if (secConstant.bending)
    {
        // 调用板单元计算
        m_plate.SetCoordinate (coord, secConstant);  //mitc3 do nothing
        Matrix Cb;
        Matrix Cs;

        // 弯曲部分
        Matrix Bb;
        Real detj;

        Matrix Bs;
        const auto nGauss = m_integForm.GetPointSize();
        for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
        {
            // secConstant.GetCb (thickness[iGauss], _OUT Cb);
            const auto& Cb6 = secConstant.m_D;
            Matrix Cb (3,3, {
                Cb6(0,0), Cb6(1,0), Cb6(2,0),
                Cb6(1,0), Cb6(3,0), Cb6(4,0),
                Cb6(2,0), Cb6(4,0), Cb6(5,0)});
            auto IOvr12 = secConstant.FI2 /12.0;
            Cb = Cb * (thickness[iGauss] * thickness[iGauss]* thickness[iGauss] * IOvr12);

            m_plate.ComputeB(iGauss, coord, _OUT Bs,_OUT Bb, _OUT detj);
            
            const Real tmp = detj * m_integForm.GetPoint(iGauss).m_weight; 
            kp += Bb.Transpose() * Cb * Bb * tmp;
            
            // secConstant.GetCs (thickness[iGauss], _OUT Cs);
            if (0 != secConstant.m_Cs.GetColSize())
            {
                auto factor = secConstant.m_ts * thickness[iGauss];
                Cs = Matrix (2, 2, {
                    secConstant.m_Cs(0, 0) * factor, secConstant.m_Cs(1, 0) * factor,
                    secConstant.m_Cs(1, 0) * factor, secConstant.m_Cs(2, 0) * factor
                });
                kp += Bs.Transpose() * Cs * Bs * tmp;
            }
            else
            {
                SGError << "Tri shell not has mid3, ";
            }

            if(secConstant.coupled)
            {
                // Matrix Cd;
                // secConstant.GetCd (thickness[iGauss], _OUT Cd); // 包含`厚度`
                const auto& B6 = secConstant.m_B;
                Matrix Cd (3,3, {
                B6(0,0), B6(1,0), B6(2,0),
                B6(1,0), B6(3,0), B6(4,0),
                B6(2,0), B6(4,0), B6(5,0)});
                Cd = Cd * (thickness[iGauss]*thickness[iGauss]);
                m_mem.ComputeStrainMatrix (coord, _OUT Bm);
                membPlateKe  += Bm.Transpose() * Cd * Bb * tmp;
            }
        }

        static const std::vector<std::size_t> indexb {2, 3, 4, 8, 9, 10, 14, 15, 16};
        keLocal.BlockFill (indexb, indexb, kp);
        
        //组装拉弯耦合项
        if (secConstant.coupled)
        {
            keLocal.BlockFill (indexm, indexb,membPlateKe);
            keLocal.BlockFill (indexb, indexm,membPlateKe.Transpose ());
        }
    }

    // 中面偏置处理
    if (fabs(zoff) > 1.0e-6)
    {
        auto TOff = SG::Algebra::identity (6);
        // const auto& zoff = cquad4.m_zoffs;
        TOff (0, 4) = zoff;
        TOff (1, 3) = -zoff;
        Matrix TToff(18, 18);
        TToff.BlockFill (0, 0, TOff);
        TToff.BlockFill (6, 6, TOff);
        TToff.BlockFill (12, 12, TOff);
        keLocal = TToff.Transpose() * keLocal * TToff;
    }
    
    kg = Tele * keLocal * Tele.Transpose();
}

/** @brief 将单元 刚度矩阵或者应力刚度矩阵从`单元局部坐标系`变换到`节点位移输出坐标系
 * `
 * @param T            3*3矩阵，   单元局部坐标系
 * @param points       单元节点在单元局部坐标系下的坐标值
 * @param dispCoordIds 单元的节点位移输出坐标系
 */
auto triCalcTransMatrix = [] (const SG::Algebra::Matrix& T,
                           const std::vector<Point>&              points,
                           const SG::DataStructure::Common::IdArray& dispCoordIds) -> Matrixd {
    Matrix Tele (18, 18);
    auto&  coordHandler = CoordinateProcessor::GetInstance ();
    for (std::size_t iNode = 0; iNode < 3; ++iNode)
    {
        const auto& cId = dispCoordIds[iNode];
        if (0 == cId)
        {
            Tele.BlockFill (iNode * 6, iNode * 6, T);
            Tele.BlockFill (iNode * 6 + 3, iNode * 6 + 3, T);
        }
        else
        {
            Matrix TDisp (3, 3);
            coordHandler.GetCoordTrans (cId, points[iNode], _OUT TDisp);
            auto tmpT = TDisp * T;
            Tele.BlockFill (iNode * 6, iNode * 6, tmpT);
            Tele.BlockFill (iNode * 6 + 3, iNode * 6 + 3, tmpT);
        }
    }
    return Tele;
};

/**
 * @brief   : 单元节点坐标，节点全局自由度索引，节点位移输出坐标系
 * @param   : nodeRepository 节点数据库
 * @param   : eleData          单元数据
 * @param   : points           节点坐标
 * @param   : globalIndex 节点位移自由度全局编号
 * @param   : dispCoordIds 位移输出坐标系变换
 * @author  : wanghuan
 * @date    : 2024.10.30
 */
void QuryElementInfo (const std::shared_ptr<SG::DBManager::INodeService> nodeRepository,
                      const CTRIA3Data&                                       eleData,
                      PointArray&                                             points,
                      std::vector<int>&                                       globalIndex,
                      std::vector<Id>&                                        dispCoordIds)
{
    auto& assembly     = AssemblyEigen::GetInstance ();
    auto& coordHandler = CoordinateProcessor::GetInstance ();

    std::vector<std::size_t> nodeIndexes (3);
    NodeData                 node;
    assembly.FindNodeIndex (eleData.m_g, _OUT nodeIndexes);

    // Node coordinate
    for (std::size_t iNode = 0; iNode < 3; ++iNode)
    {
        nodeRepository->FindByIndex (nodeIndexes[iNode], _OUT node);
        coordHandler.GetGlobalCoordinate (node, _OUT points[iNode], _OUT dispCoordIds[iNode]);
    }

    static const DofSet DOFs{ DOFDirections::DOF_Ux,   DOFDirections::DOF_Uy,   DOFDirections::DOF_Uz,
                              DOFDirections::DOF_Rotx, DOFDirections::DOF_Roty, DOFDirections::DOF_Rotz };
    assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT globalIndex);
}

void Tria3ShellCalculator::ComputeMass
(   const Real&                area,
    const ShellSectionProp1&                       secConstant,
    //偏置计算存在问题，暂设为可修改左值
    SG::DataStructure::Common::Real                 zoffs,    
    _OUT SG::Algebra::Matrixd&                coupledM,
    _OUT SG::Algebra::Matrixd&                lumpedM
)
{   
    Real RHOt = secConstant.m_density *secConstant.m_t;
    Real m = (RHOt + secConstant.m_nsm)* area;

    //强制修改：目前偏置计算不正确先令zoffs = 0.0，以通过无偏置情况的计算
    zoffs = 0.0;
    //偏置质心发生变化带来的修正值
    Real m_c = m * zoffs * zoffs;

    // Real I_x = (m / 36.0) * ((localCoord(0, 1) - localCoord(1, 1)) * (localCoord(0, 1) - localCoord(1, 1)) + (localCoord(1, 1) - localCoord(2, 1)) * (localCoord(1, 1) - localCoord(2, 1))
            //  + (localCoord(0, 1) - localCoord(2, 1)) * (localCoord(0, 1) * localCoord(2, 1)));
    // Real I_y = (m / 36.0) * ((localCoord(1, 0) - localCoord(2, 0)) * (localCoord(1, 0) - localCoord(2, 0)) + (localCoord(2, 0) - localCoord(0, 0)) * (localCoord(2, 0) - localCoord(0, 0))
            //  + (localCoord(0, 0) - localCoord(1, 0)) * (localCoord(0, 0) - localCoord(1, 0)));
    
    Matrix templumpedM(6, 6);
    templumpedM(0, 0) = m / 3.0 + m_c;
    templumpedM(1, 1) = m / 3.0 + m_c;
    templumpedM(2, 2) = m / 3.0 + m_c;
    // templumpedM(3, 3) = I_x /3.0;
    // templumpedM(4, 4) = I_y /3.0;
    // templumpedM(5, 5) = .0;
    lumpedM.BlockFill (0, 0, templumpedM);
    lumpedM.BlockFill (6, 6, templumpedM);
    lumpedM.BlockFill (12, 12, templumpedM);

    //偏置对耦合项带来的修正值
    Matrix e1 = SG::Algebra::identity (3) * (m / 12.0) + SG::Algebra::identity (3) * (m / 12.0) * zoffs * zoffs;
    //偏置质心发生变化带来的修正值
    Matrix e2 = SG::Algebra::identity (3) * (m / 12.0) * 2.0 +SG::Algebra::identity (3) * (m / 12.0) * zoffs * zoffs * 2.0;
    //偏置对转动项带来的修正值
    Matrix r1 = SG::Algebra::identity (3) * (m / 12.0) * zoffs * zoffs;
    Matrix r2 = SG::Algebra::identity (3) * (m / 12.0) * zoffs * zoffs * 2.0;
    Matrix tempcoupledM1(6, 6);
    Matrix tempcoupledM2(6, 6);
    tempcoupledM2.BlockFill (0, 0, e2);
    tempcoupledM2.BlockFill (3, 3, r2);
    tempcoupledM1.BlockFill (0, 0, e1);
    tempcoupledM1.BlockFill (3, 3, r1);
    coupledM.BlockFill (0, 0, tempcoupledM2);
    coupledM.BlockFill (6, 6, tempcoupledM2);
    coupledM.BlockFill (12, 12, tempcoupledM2);
    coupledM.BlockFill (6, 0, tempcoupledM1);
    coupledM.BlockFill (12, 0, tempcoupledM1);
    coupledM.BlockFill (0, 6,  tempcoupledM1);
    coupledM.BlockFill (0, 12, tempcoupledM1);
    coupledM.BlockFill (6, 12, tempcoupledM1);
    coupledM.BlockFill (12, 6, tempcoupledM1);
}

/** @brief 板单元计算并保存刚度矩阵
 * 
 */
void Tria3ShellCalculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto& assembly = AssemblyEigen::GetInstance();
    const auto& nodeRepository = dbServiceFactory -> get<INodeService>();
    const auto& eleRepository = dbServiceFactory -> get<IElementService>();

    // 壳属性计算器
    ShellSectionParser1 parser;
    ShellSectionProp1 secConstant;
    CTRIA3Data eleData;
    std::vector<Point> points(3);
    std::vector<int>   globalIndex (18);
    std::vector<Id> dispCoordIds (3);
    Matrix kGloal (18, 18);
    Matrix LumpedM(18, 18);
    Matrix CoupledM(18, 18);

    static const DofSet DOFs{ DOFDirections::DOF_Ux, DOFDirections::DOF_Uy, DOFDirections::DOF_Uz };
    
    const auto eleSize = eleRepository -> Size<SG::DataStructure::FEM::CTRIA3Data>();
    for (std::size_t index{0}; index < eleSize; ++index)
    {
        /** 从数据库获取对应的单元数据 */
        eleRepository -> FindByIndex(index, _OUT eleData);
        QuryElementInfo (nodeRepository, eleData, _OUT points, _OUT globalIndex, _OUT dispCoordIds);

        // Local Coordinate system
        Matrix T, localCoord;
        Point origin;
        ElementTools::computePlaneCoordTrans (points, _OUT T, origin, _OUT localCoord);
        
        
        // Section Property
        Status status = parser.Compute (dbServiceFactory, eleData.m_pId, _OUT secConstant);
        assert (Status::STATUS_Success == status);
        if(secConstant.needTrans)
        {
            // 单元形心点
            Point eleCentral{ (points[0].m_x + points[1].m_x + points[2].m_x) /3.0,
                (points[0].m_y + points[1].m_y + points[2].m_y) /3.0,
                (points[0].m_z + points[1].m_z + points[2].m_z) /3.0 };
                
            // G1-G2节点向量:
            vector<Real> vectorG1G2 = { localCoord (1, 0) - localCoord (0, 0), localCoord (1, 1) - localCoord (0, 1) };
            assert (fabs(vectorG1G2[1]) < 1.0e-6 ); // 单元局部坐标系x轴与 1-2节点连成的边重合
                
            // 单元材料铺层方向
            vector<Real> eleMatAxis;
            Real         theta      = eleData.m_theta * PI / 180.0;
            bool         transFlag  = SG::FEMSolver::Element::calcMaterialOritation (eleData.m_mcId, eleCentral, T, vectorG1G2, theta, _OUT eleMatAxis);
            if (transFlag)
            {
                secConstant = secConstant.Translate (eleMatAxis[0], eleMatAxis[1]);
            }

        }

        // 刚度矩阵变换到全局坐标系
        auto Tele = triCalcTransMatrix (T, points, dispCoordIds);

        const auto zoff = eleData.m_zoffs + (secConstant.m_z0 + 0.5* secConstant.m_t);
        ComputeStiffness (localCoord, secConstant, zoff, Tele, _OUT kGloal);
        assembly.AccumulateStiffness (globalIndex, kGloal);

        // 三角形面积计算
        //A = (x1 * (y2 - y3) + y1 * (y3 - y1) + x2 * (y1 - y2)) / 2.0;
        Real A = (localCoord(0, 0) * (localCoord(1, 1) - localCoord(2, 1)) + localCoord(1, 0) * (localCoord(2, 1) - localCoord(0, 1)) + localCoord(2, 0) * (localCoord(0, 1) - localCoord(1, 1))) / 2.0;

        ComputeMass (A, secConstant, zoff, _OUT CoupledM, _OUT LumpedM);
        const auto coupled = assembly.NeedCoupleMass ();
        if (coupled)
        {
            CoupledM = Tele.Transpose() * CoupledM * Tele;
            assembly.AccumulateCoupledMass (globalIndex, CoupledM);
        }
        else
        {
            assembly.AccumulateCoupledMass (globalIndex, LumpedM);
        }
    }
}

void Tria3ShellCalculator::DataRecovery(const SG::DBManager::DBServiceFactorySP dbServiceFactory, int subcaseId, SG::DataStructure::Common::Real time, const std::vector<SG::DataStructure::Common::Real>& globalDisplacement) 
{
    // SG::FEMSolver::Element::Gauss2DTria3P m_integForm;  ///< 积分格式（三点积分）

    // //单元节点全局坐标系位移向量
    // std::vector<SG::DataStructure::Common::Real> eleDisplacement(18);
    // //单元节点单元局部坐标系位移向量
    // std::vector<SG::DataStructure::Common::Real> eleLocDisplacement(18);
    // //为匹配膜、板弯曲、板剪切应变矩阵维度，将壳位移拆分成膜位移及板弯曲位移
    // //单元节点膜单元局部坐标系位移向量
    // std::vector<SG::DataStructure::Common::Real> eleMemDisplacement(9);
    // //单元节点板单元局部坐标系位移向量
    // std::vector<SG::DataStructure::Common::Real> elePlateDisplacement(9);

    // CTRIA3Data eleData;
    // // std::vector<std::size_t> nodeIndexes(3);
    // std::vector<int>dofGlobalIds(18);
    // // NodeData node;
    // std::vector<Point> points(3);
    // std::vector<Id> dispCoordIds (3);
    // //工况信息
    // SG::DataStructure::FEM::SubCase subcase;
    // // SG::FEMSolver::Element::Plane3NodeShpFunc m_shpFunc; ///< 积分点上的形函数

    // // 壳属性计算器
    // ShellSectionParser parser;
    // ShellSecMaterialConstant secConstant;


    // //获取后处理输出变量
    // dbServiceFactory->get<SubCaseService>()->FindById(subcaseId, subcase);
    // const auto& nodeRepository = dbServiceFactory -> get<NodeService>();
    // // OutputSetting outPutSetList=subcase.m_out;

    // //输出积分点应力
    // if(true)
    // {
    //     //积分点应力存储Vector
    //     std::vector<SG::Algebra::Matrixd> stressSet;

    //     const auto eleSize = eleRepository -> Size();
    //     for (std::size_t index{ 0 }; index < eleSize; ++index)
    //     {
    //         // 单元积分点在单元坐标系下的应变矩阵行为分量依次为：xx,yy,zz,xy,yz,xz，列为积分点
    //         Matrix eleIntegStrain (6, 3);
    //         // 单元积分点在单元坐标系下的应力矩阵行为分量依次为：xx,yy,zz,xy,yz,xz，列为积分点
    //         Matrix eleIntegStress (6, 3);
    //         // 根据索引遍历读取单元信息
    //         dbServiceFactory -> get<CTRIA3Service>()->FindByIndex(index, _OUT eleData);
    //         QuryElementInfo (nodeRepository, eleData, _OUT points, _OUT dofGlobalIds, _OUT dispCoordIds);

            
    //         //将单元全局坐标系下的节点信息转换到单元局部坐标系，存储于localCoord
    //         Matrix T, localCoord;
    //         computePlaneCoordTrans (points,_OUT T, _OUT localCoord); 
            
    //         //将单元节点位移从全局坐标系转换到单元局部坐标系
    //         auto Tele = triCalcTransMatrix (T, points, dispCoordIds);
    //         Tele = Tele.Transpose ();

    //         //从位移向量中读取对应节点位移
    //         for(int idof=0;idof<18;idof++){
    //             eleDisplacement[idof]=globalDisplacement[dofGlobalIds[idof]];
    //         }
        
            
    //         // 全局坐标系位移转换为单元坐标系位移
    //         //未重载矩阵*向量，此处先用for循环替代
    //         for(int i=0;i<18;i++)
    //         for (int j=0;j<18;j++)
    //             eleLocDisplacement[i]+=Tele[i][j]*eleDisplacement[j];

    //         //单元节点膜单元局部坐标系位移向量
    //         std::vector<SG::DataStructure::Common::Real> eleMemDisplacement(9);
    //         //单元节点板单元局部坐标系位移向量
    //         std::vector<SG::DataStructure::Common::Real> elePlateDisplacement(9);
    //         //拆分单元位移
    //         for(int i=0;i<3;i++){
    //             eleMemDisplacement[3*i]=eleLocDisplacement[6*i];
    //             eleMemDisplacement[3*i+1]=eleLocDisplacement[6*i+1];
    //             eleMemDisplacement[3*i+2]=eleLocDisplacement[6*i+5];
    //             elePlateDisplacement[3*i]=eleLocDisplacement[6*i+2];
    //             elePlateDisplacement[3*i+1]=eleLocDisplacement[6*i+3];
    //             elePlateDisplacement[3*i+2]=eleLocDisplacement[6*i+4];
    //         }

    //         // 膜应变矩阵Bm，板弯曲应变矩阵Bb，板剪切应变矩阵Bs
    //         Matrix Bm, Bb, Bs;

    //         const auto nGauss = m_integForm.GetPointSize();
    //         for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss){

    //             Real jm, jb;
    //             // 调用Membrane单元进行计算协调部分应变矩阵Bm
    //             m_mem.ComputeStrainMatrix (iGauss, localCoord, Bm, jm);
    //             //计算膜部分应变：
    //             for(int i=0;i<4;i++){
    //             eleIntegStrain[0][iGauss]+=Bm[0][3*i]*eleMemDisplacement[3*i];
    //             eleIntegStrain[1][iGauss]+=Bm[1][3*i+1]*eleMemDisplacement[3*i+1];
    //             eleIntegStrain[3][iGauss]+=(Bm[2][3*i+1]*eleMemDisplacement[3*i+1]+Bm[2][3*i]*eleMemDisplacement[3*i]);
    //             }
    //             // 调用板单元计算
    //             m_plate.ComputeBb (iGauss, localCoord, _OUT Bb, _OUT jb);
    //             m_plate.ComputeBs (iGauss, jb, localCoord, _OUT Bs);
    //             //计算板剪切应变（根据屈曲计算需求，此处计算中面）
    //             for(int i=0;i<4;i++){
    //                 // 应变矩阵Bs顺序与常用排序相反，在此调整
    //                 eleIntegStrain[5][iGauss] += Bs[0][3 * i] * elePlateDisplacement[3 * i]
    //                                                 + Bs[0][3 * i + 1] * elePlateDisplacement[3 * i + 1]
    //                                                 + Bs[0][3 * i + 2] * elePlateDisplacement[3 * i + 2];
    //                 eleIntegStrain[4][iGauss] += Bs[1][3 * i] * elePlateDisplacement[3 * i]
    //                                                 + Bs[1][3 * i + 1] * elePlateDisplacement[3 * i + 1]
    //                                                 + Bs[1][3 * i + 2] * elePlateDisplacement[3 * i + 2];
    //             }

    //         }
    //         Matrix elasticMatrix (6, 6);
    //         for (int i = 0; i < 2; i++)
    //             for (int j = 0; j < 2; j++)
    //                 elasticMatrix[i][j] = secConstant.m_Cm[i][j] / secConstant.m_t;
    //         elasticMatrix[3][3] = secConstant.m_Cm[2][2] / secConstant.m_t;
    //         elasticMatrix[4][4] = secConstant.m_Cm[2][2] * 5.0 / 6.0;
    //         elasticMatrix[5][5] = secConstant.m_Cm[2][2] * 5.0 / 6.0;
    //         // 存储单元局部坐标系下积分点应力，格式与积分点应变相同
    //         eleIntegStress = elasticMatrix * eleIntegStrain;
    //         stressSet.push_back (eleIntegStress);
    //         // for(int i =0;i<3 ;i++)
    //         // {
    //         //     for(int j=0;j<6;j++)
    //         //     {
    //         //         std :: cout << eleIntegStrain[i][j]<< "  " ;
    //         //     }
    //         //     std :: cout <<std ::endl;
    //         // }
    //     }

    //     dbServiceFactory->get<CTRIA3IntegStressService>()->add(subcaseId, stressSet);
    // }
}

void Tria3ShellCalculator::ComputeStressStiffness(const SG::DBManager::DBServiceFactorySP dbServiceFactory,int subcaseId)
{

}